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r 

HETEROGENEOUS MEDIA FOR NEUTRON ‘DIFFUSION CALCULATIONS 


A. Jacobs 


I. INTRODUCTION 


We shall develop, in this work, the mathematical formulation 
of a new method of treating neutron diffusion in certain types of 
heterogeneous media. Hie heterogeneity of immediate concern, and 
toward which this work is directed, is that of a regular array of 
vacuum channels (such as a square lattice of cylindrical holes) in 



k 


» 



an otherwise homogeneous medium. With modification of details, 
the general procedure should be applicable to other types of heter 
ogeneity. However, a requirement which should be imposed is that 
the heterogeneity results in two characteristic directions. For 
example, in the case of a regular array of vacuum channels, the 
two characteristic directions are parallel and perpendicular to 
the channel axis. 

Probably the most important considerations of neutron dif- 
fusion In media with holes are those of Behrens (B). However, all 
previous work, including (B), is devoted to the determination of 
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the effects, of heterogeneity on specific parameters relevant to 
neutron diffusion, rather than to a formulation of a general method 
from which various descriptions of the neutron distribution can be 
determined- The general problem may be stated along these lines: 

In media with heterogeneity, such as the type considered here, it 
is clear that a "simple homogenization" of the medium for neutron 
diffusion calculations is not a valid representation. We define a 
simple homogeni ration as the process of reducing the medium cross 
sections merely by the ratio of material volume to material-plus- 
vacuum volume. The streaming of neutrons in the vacuum channels 
leads to a spreading of the neutrons in the longitudinal direction 
(parallel to channel axis) which is larger than that in the trans- 
verse direction (perpendicular to channel axis). Thus, we find 
that not only is the simple homogenization questionable for omni- 
directional parameter calculation, but also the anisotropic effects 
are completely subdued. Of course, the reason that we have for con- 
sidering homogenization of the medium is the existence of an "arsenal" 
of possible mathematical, attacks for such problems. 

The new approach to homogenization, which we shall develop, 
is based on the reasoning that neutrons traveling with a large 
component of their velocity in the longitudinal direction probably 
travel further between collisions, on the average, than those travel- 
ing with a large component of velocity in the transverse direction, 
let us introduce this effect into the neutron transport equation for 


the homogenized medium by allowing the mean free path to be angular- 
dependent. Clearly, by bo doing, we have embarked on a mathematical 
fiction since the mean free path, as used in neutron transport cal- 
culations, is a local parameter. This concept is certainly no more 
confusing than the idea of medium homogenization. Let us stress 


that we are Imposing a total cross section which varies according 
to the direction of neutron travel and not merely the usual scat- 
tering directional dependence. 

All considerations will be based on the idealization of 
monoenergetic neutron transport. Although time -dependent equations 
will be developed, most of the analysis will be devoted to the 
calculation of stationary states. For the purpose of completeness, 
we shall develop the mathematical formulation along several lines. 
Sp^ifically : We shall consider both the neutron flux and collison 

density as dependent variables. We shall apply the familiar and 
double-P^ approximations, as well as the moment decomposition. We 
shall demonstrate that, for the case of isotropic scattering, the 
normal mode procedure, recently used for the solution of several 
types of neutron transport problems, is applicable and yields exact, 
closed -form solutions. 

Ihe major part of this work we direct toward the mathematical 
formulation and physical interpretation of the new theory. We shall, 
however, devote a section to brief remarks relevant to applica- 


tion. There is meager experimental data available for lattices of 


the type considered. It will therefore be difficult to evaluate 
the present theory. These considerations will yield a well-defined, 
albeit not veil-substantiated, route to solution of problems involv- 
ing neutron diffusion in media pierced by vacuum channels. 

II. MATHEMATICAL FORMULATION 

It would seem that application of these ideas to finite media 

dictates the use of non-separable position-angle -dependent mean free 

% 

paths. The inclusion of position-dependence leads to gross diffi- 
culties which we have as yet not resolved. We shall assume that 
the mean free path in the homogenized medium depends only on the 
angle between the neutron velocity and the direction of the position 
variable. In all calculations ve assume plane symmetry with the 
position variable either along the longitudinal or transverse 
direction. 

2. 1 The Neutron Flux Equation 

The monoenergetic neutron transport equation for homogeneous 
media with plane symmetry may be written as 

| + (jt \|r(x,(i,t) + o(u) *(x,u,t) 

* C J a(n») f(g*2») *(x,p',t) dfl* + S(x,p,t) 

where y(x,p,t) is the neutron flux distribution as a function of 
position, x, direction cosine of neutron travel relative to 
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x-direction, p, and time, t; v 1 b neutron speed; 0(4) 1 b the 

■y 

angular -dependent total cross section; c is the mean number of 
secondary neutrons which emanate from a neutron -nucleus collision; 

is the normalized distribution in fl, the neutron post-collision 
direction of travel, of secondary neutrons produced by collision of 
a primary neutron with pre-collision direction and, S(x,u,t) is 
the rate of neutron introduction from sources which are independent 
of the neutron distribution. Although c and f are not necessarily 
descriptive < of a non-multiplying medium, we shall use the terms 
scattering probability for c and scattering distribution for f to 
avoid s^il^ed discourse. We employ an expansion of the scattering 
distribution in terms of Legendre polynomials |p n (£*fi , )| , i.e.. 


r(n-O') 



2n + l, 
~TT~ f n 


p n (£?'0 r ) 


(2.1.2) 


and the spherical harmonics addition theorem (H, p. 1^3) to eliminate 
the azimuthal direction dependence appearing in the integral in 
Eq.. (2.1.1). WO obtain 

~ t(*>4>t) + 4 tU,u,t) + cr(u) *(x,u,t) 

♦ 1 

= I £ < 2n + f n f P n (y,) a ^' } t( x ^',t) V 

n •'-1 


+ S(x,n,t) 


(2.1.3) 
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* 


In order to proceed with a general discussion of the 
properties of Eq. (2.1.3) we require a more specific representation 
of 0 ( 4 ). We suppose, with little loss in relevant generality, that 

* £ °n P „<“> 
n 

Let us further define the sets j^ n (x,t)| and |S n (x,t)| by 

* 1 

♦ n (x>t) = f p (n) ♦(x,u, t) dp 
•'-1 

*1 

S Q (x,t) = J P n ( u ) S(x,p,t) dn 

- 1 

Equations (2.1.5a) and (2.1.5b) specify the respective expansion 
coefficients in Legendre polynomial expansions of the neutron flux 
and source density, e.g. , 

+(x,u,t) = 2 — * n (x,t) P (p) 

n 


In these terms, integration of Eq. (2.1. 3 ) over the interval 
p€ ( —1, +l) yields the relation 


“5t + 


^1 

~5x + 


(1 



a dr 
n T n 


S 0 


Eq. (2.1.7) is the continuity equation for neutron motion. The 
only term which appears in an unfamiliar fbrm is that which expresses 
the total neutron interaction rate. Clearly, 


♦ 1 

Y. a n ^n( x,t ^ -/ 0 (u) t(*,u,t) d|i 


(2.1.4) 

(2.1.5a) 

(2.1.5b) 

( 2 . 1 . 6 ) 

(2.1.7) 


( 2 . 1 . 8 ) 
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A further, familiar, reduction of the transport equation 
can be made in terms of the sets {S q J, and | • Using the 

recurrence relation for Legendre polynomials (H, p. 32), 


(2n + 1) n P n ( u ) = n (u) + (n + 1) P^ (u) 


(2.1.9) 


yields the set of coupled differential equations. 


1 d 

v 57 




n 


2 n 




n 4 


1 5 


2n 4 1 


5x ^n4i^ x '^ 


+ {2t 4 X) °m A 4nn (l " cf n } 

t, m 

« S n (x,t) (2.1.10) 


Hie set 



is defined by 



P z (») P ( M ) P n (u) dp 


( 2 . 1 . 11 ) 


and has the following properties (H, p. 87 ): 

(i) Hie order of the indices is unimportant. 

(ii) A Zrm = 0 if the sum of any two of the indices 

is less than the third. 

(iii) A f mr) = 0 if ( 4 b 4 n is odd. 

(iv) When A^ mn £ 0, i.e., avoiding (ii) and (iii), we have 
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Vmn + m + n + lj 


+ m -.el- 1 ! 

(2)(4)---(4 + m - n) 


0 


(l)f3) — C-t + n - a - 1) (1)(3) — (n + a - l - 1) 

(2)t4).:.(i - n - a) 'J Ln2)TO-'-Tn“»'"' ; Zr^. 


A2)W" ■{£ + a + n) 

( 1 ) (3 ) * ••(■£ + m + n - 


iT 

It is from Eq. (2.1.10) that the various approximations to 
the transport equation are derived. Eefore we proceed with detailed 
examination of some approximations, let us reexamine Eq. (2.1.3). 

We obtain a simpler integral term if the dependent variable is 
changed to the neutron collision density, 

1 

F(x,u,t) = a(n) t(x,u,t) . 

2.2 The Collision Density Equation 

We reformulate Eq. (2.1.3) in terms of the collision 
density, F(x,p,t), and the mean free path, A(p) = l/o(p), and 
obtain 

F (|c,u,t) + p A((i) ^ F(x,p,t) + F(x,u,t) 




I X) ( 2n + f n P n (M) f P n^') F (x,n*,t) dp' 

n J-& 


. 1 . 12 ) 


+ S(x,u, t ) 


( 2 . 2 . 1 ) 
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With the Bets |F Q (x,t)| and \^ n \ defined by 


F(x,»i,t) - ^_Li F Q (x,t) P Q (u) 


(2.2.2a) 


Mu) “ S X n P n (u) 


(2.2.2b) 


we obtain the set of coupled differential equations (cf., Eq. (2. 1. 10) ) , 

S in * 1) T 5t V x,t> + £ (a + h V x>t) 

+ 1TTT < a + 57 V*’* 5 


+ (1 - c f n ) F n (x,t) 


= S Q (x,t) 


(2.2.3) 


We note that if A =0 for n > 0, i.e., the familiar case of an 
n 

angle -independent mean free path, then Eq. .(2.2.3), with the aid of 
the properties of j, reduces to 

N>». . nX 0 S *n-1 + (o * \ 8f n.l . . , 

v 5t“ 2n + 1 dt 2n + 1 dt n □ 


which is the expected result. We note further that the n 
member of Eq. (2.2.3) yields the continuity equation 

_ dF OF 

2-r \ ZT + 2-r (2t + U \ A £m i ^ + (i - c ) F Q ■= S Q 
m 


- o 


(2.2.4) 


(2.2.50 


S 
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In recognizing Eq. (2.2.5) as the continuity equation, we have 
used the easily derived relations 

- E K V 1 ** 1 

n 

y, (x,t) = Y, (2n + i} \ A nml F n (x ' t} 
m, n 

In the remainder of this work we shall assume the case of a 
stationary state. In most cases this assumption merely leads to 
simpler algebra and notation and is actually not a requirement for 
the determination of a solution. We shall discuss the P^-approxima- 
tion and double -P -approximation as applied to both flux and 
collision density expansions. We shall consider the moment decom- 
position for both flux and collision density. And, finally, we 
shall consider, in greater depth, the case of isotropic scattering 

using the collision density as dependent variable. 

1 

2.3 The P^-Approximation 

We define the P N ~approximation based on a flux expansion, 
or collision density expansion, by the requirements that in 
Eq. (2.1.10), or Eq. (2.2.3)/ ^ n {x) = ® Tor o > N, or F n (x) = 0 f.or 
n > N, and the equations labelled by n > N are discarded. Thus, 

In a P^-approximation we obtain H + 1 coupled differential equations 

» 

with the N + 1 dependent variables ^ n (x), n = 0, 1, • • • , N, or 

t 


(2.2.6a) 

(2.2.6b) 
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F r (x), n = 0, 1, •••, N. For example, the P^-approxiaation based 
on a flux expansion gives the single relation 


°0 V 1 “ c ) *q(*) = 0 

which has the implication c = 1 for non-trivial to * KrLs is a 
familiar result* The ? D -approximation based on a collision density 
expansion gives the unusual relation 


X 1 ^0 
3 dx + (' 


c)F q (x) = 0 


After this present brief comment we shall restrict our consideration 
to the case of A(p) a symmetric function of 4 on the interval 
pe(-l, +1), and, in that case Ai = 0. For the moment, let us sup- 
pose that \ f 0. To be specific, we suppose that A^ > 0. 

According to Eq. (2.3*2), the implication is that cLF^dx < 0 which 
indicates a neutron flow in the + x - direction. The question which 
we must pose iss Can F be 4 -lndependent (consistent with P o -approxi- 
raation) and yet have ^(x,p) represent a neutron flow in + x - direction 
(i.e., ^ increase with 4 )? Clearly the answer is in the affirmative 
since, if A(p) increases with u (implied by Aj > 0 ), then the ratio 
\Jr/ A can be 4 -independent if f also increases with 4. In all follow- 
ing considerations we assume that A( 4 ), and thus 0 ( 4 ), is a 
symmetric function of 4 . Thus we shall always require A^ = 0, and 
a n = 0 , for n odd. 


(2.3-1) 


(2.3-2) 


- 12 - 

The Pj -approximation (i.e., diffusion theory) based on a 
flux expansion gives the two equations 

d t, 

S” + (1 - c) o 0 t 0 (x) » S 0 (x) 
d^r 

J * (! - cf i)( a o + 5 a 2 ) ♦i(x) = S 1 (x) 

If it is further assumed that all neutron sources are isotropic 
such that S q (x) * 0 for n > 0, then Eqs. (2-3*3 a ) and (2.3*3^) com- 
bine to give the usual diffusion theory relations 

1 0 

- D — r— + 0* t„(x) « S (x) 
dx 

, s d *o 

'iW * - D dsr 

with diffusion coefficient, D, and "absorption" cross section, o', 
.given by 

D = 3(1 - cf 1 )(o 0 + 2 c 2 /5) 

O' = (1 - C) (Jq 

1 

It 6hould be noted that only o 0 and a 2 enter into the diffusion 
theory parameters. The o(u) expansion was not truncated at a 
quadratic in order to obtain Eqs. (2. 3- 5a) and (2. 3. 5b). In passing, 
it is also worth noting that, had we included time -dependence, the 


(2.3.3a) 

(2.3.3b) 

(2.3.4a) 

(2.3.4b) 

( 2 . 3 . 5 a) 

( 2 . 3 . 5 b) 


I 


-13- 


Pl -approximation would have given the "telegraphist's equation." 

Ihe added assumption that ctyj/ v St < < S^q/Sx results in the form 
of the familiar time -dependent diffusion theory with D and a' again 
given by Eqs. (2. 3- 5a) and (2.3.5b). 

Hhe Pj -approximation based on a collision density expansion 
gives the two equations 

2 dF l 

(\ + j \ ) -g- + (l-o) Fq (x) = S D (x) 

1 2 ^0 

+ f V 5T + ^ - cf i> F 1^) - SiW 

In the case of isotropic sources, Eqs. (2.3>6a) and (2.3.6b) combine 
to give 

d 2 F 

_ D . _ + (1 . c ) F (x) = S (x) 
dx 

. dF 

F (x) = — 7 ^ t — 

1 V } (\ + 2 \/ 5) dx 

where the "diffusion coefficient" is 

(\ + 2N/5) > 

~ 3(1 - cfj) 


(2.3.6a) 

(2.3.6b) 

(2.3. T«) 
(2.3.7b) 

(2.3.8) 


As is the case in the flux based diffusion theory, a quadratic 
truncation of A(p) is not necessary to obtain the results given in 
Eqs, (2.3.7) and (2.3.8). 
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The usual result of exponential spatial decline away from 
sources in infinite media is found for flux and collision density. 
Hie "period" of the exponential, the diffusion length L, is given 
by 

j2 1 

* * 3o 0 (l^ c)(l - cfj)(o 0 + 2a 2 /5) 


based on a flux expansion, and 



(\ * 2>> 2 /5) 2 

3(1 - c)!l - cf J J 


(2.3.9) 


(2.3.10) 


based on a collision density expansion. We have pointed out the 
dissimilar results obtained from a P 0 -approximation based on flux 
and collision density expansions. We can illustrate the divergence 
of results for the Pj -approximation by considering the ratio L^/Lp 
for a given problem. Let us suppose that A(p) is actually an even 
quadratic, i.e., = 0 for n = 1 and n > 2. The corresponding 0 ( 4 ) 

is not a quadratic, but we need compute only o 0 and a 2 since these 
coefficients determine L^.. In Figure 1 we display the results for 
the range X 2 /A q €( 0,2). Clearly, the flux and collision density based 
expansions can lead to significantly different results. 

It should be expected that the two diffusion theories give 
different results. We note that whereas the ^-approximation yields 
an accurate representation for neutron current, ijr^x), but an inaccurate 
(truncated) representation for collision density, the F-approximation 
results in an accurate representation for total interaction rate. 



I 


f 
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F^(x), but an inaccurate representation for current. For several 
reasons it will seem that we favor the collision density expansion 
in this work. However, these reasons are mainly of an algebraic 
character and it should not be construed that the F-formulation is, 
in all cases, superior. 

Approximations of higher order than Pj are accomplished 
following the usual general prescriptions. The added complications 
due to the angular-dependence of the mean free path place no restric- 
tions on the formalism. Higher order approximations lessen the 
differences exhibited by the ijr and F-formulations. 

2.4 The Double -P^-Approxlmat ion 

The double -P^-approximation is derived from Yvon's 
method whereby the angular-dependence of the neutron flux, or col- 
lision density, is decomposed into contributions from + x - directed 
and from - x - directed neutrons. In order to simplify notation, 
let us consider the case of a stationary state in a medium charac- 
terized by isotropic scattering. The plane symmetry, monoenergetic 
neutron transport equation, Eq. (2.1.3), is then 

‘ /' 

4 t(k,u) + a(p) *(x,p) = | / o(p') tKx,u') dp* + S(x,p) 

-1 


(2.4.1 
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We use the half -angle -range expansions 


= X) (2n + 1) **(x) P Q (2p - 1), M > 0 

n 

*53 (2n + 1) *~(x) p n (2u +1), 4 < 0 

n 

S(x,u) = E (2n + 1) S*(x) P n (2 M - 1), n > 0 

n 

= 53 (2n + 1) S~(x) P n (2 M + 1),. n < 0 

n 

o(m) = 53 °n P Q ( 2 4 “ 1 )> p > 0 

n 

* 53 P n (2M + l)f M < 0 

n 


to obtain the set of coupled differential equations 

* 2 T) (21 * 1) o* 

4,m 

i *? 4 K 


n 


d\|r 


n-1 


2n + 1 dx 


d\lr d\Jr~ 

n + 1 n+1 ± __r 

2n + 1 dx dx 


= c 6 


no 


2>i 


*4 4 


Hie double-P^-approximation is defined by the requirement 
\jr*(x) = 0 for n > H and the equations labelled by n > N are 


(2. 4.2a) 
(2.4.2b) 
(2.4.2c) 
(2. 4. 2d) 

(2. 4. 2e ) 

(2.4.2f) 

(2.4.3) 


discarded. 
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Ihe same analysis is followed for the collision density 
based approximation. Thus, with the definitions 


n 


£ 

n 


(2n + 1) F*(x) P n (2p - 1), 

4 > 0 

(2.4.4a) 

(2n + 1) F"(x) P (2u + 1), 
n n 

4 < 0 

(2.4.4b) 

K P n‘^ - V ’ 

4 > 0 

(2.4.4c) 

K p „(^ + 1). 

4 < 0 

(2. 4. 4d) 


ve obtain the transport equation in the form 


2n 


TT £ (2t . 1) V A &j0 _ 


fi 

1 dx 


dF. 


n + 1 V* / . , + V4i » 

+ 2n + 1 {2C + 1) \ A m,n+1 


V* + dF „ + 

P {U * 1) X A { -£ ♦ 2F 

<-,m 


c ( F n + F ii) + 2S± 

VQ 0 0' n 


(2.4.5) 


The collision density double-P^-approximation is defined as in 
the y -formulation. 
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Besides the usual usefulness of Yvon's method in the 
solution of problems where accurate representation of source or 
free boundaries are required, we find an added flexibility for the 
angular -dependence of the mean free path. Thus, A(n) may have one 
form for 4 > 0 determined by the set |a*^ , or , and another 

form for 4 < 0 following the set | A~|, or |o~ j . For example, 
using these methods we can express a symmetric A(u) which varies 
linearly with 4 for 4 > 0 and 4 < 0 by using only two terms in each 
A(u) expansion, i.e., set A Q = A Q , Aj = — Aj , and A~ = 0 for 
n > 1 . 


2. 5 A Moment Decomposition 


In the usual theory of neutron transport through homo- 
geneous media, it is well-known that any space -angle moment of the 
neutron distribution can be found even though the distribution itself 
is unknown. In fact, an important method of determining the neutron 
distribution is the construction of a likely flux shape from a 
finite set of moments. Let us now consider a moment decomposition 
for the case of an angular-dependent mean free path where, as in our 
previous considerations, two formulations will be examined, i.e., 
with \jr and F as dependent variable. 

We define the neutron flux and source moments by 



+•* 



(x) dx 
n' ' 


S n (x) dx 


(2.5-la) 


(2.5.1b) 
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We assume that the medium is of infinite extent. Multiplication 
of the stationary state form of Eq. (2.1.10) by x^ and integration 
over x C (-», +•) yields the set of algebraic moment relations 

^ Oa - 1) % (i - cf n ) tj - sJ . [ nt J : > 4 ( „ 4 


With the collision density moments similiarily defined, i.e., 

+m 

F ^ = / x^ F (x) dx 

* n L n 


( 2 - 5 - 2 ) 


(2-5-3) 


and performing similar operations on the stationary state form of 
Eq. (2.2.3 )t( ve find a set of algebraic equations relating the 
collision density moments, i.e., 


(1 - cf ) F* 3 = S J + 6 ^ (2 1 + 1) A A. ,F^* 1 

v n 7 n n 2n + 1 ' m -cm,n-| £ 


+ \ ' ( 2 / + 1) A A» F" 3- * 

2n + 1 4 ^ ' ' m tm,n+l £ 

■c,m 


(2.5-4) 


Etie moments of the neutron distribution resulting from a unit, 
plane, isotropic source (at x = 0) are easily interpre table in terms 
of important macroscopic parameters. For this source, = 5^ 6j 0 . 
Let us first consider calculation of flux moments and then contrast 
this result with the method applied to collision density moments. 

For the sake of definiteness, let us assume that 0(4) is an even M 
degree polynomial in 4. With 0(4) an even function, it is clear that 
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^ = 0 for j + n odd. An examination of Eq. (2.5*2), using the 
properties of the set | A /mn | , indicates that with M ^ 0 and n + j 

even, ^ depends on ^ , *{ n _ Mr t*j n . M|+a >***> *n +M _, ' ^n+M * 

Thus, we find that we cannot determine without employing a trunca- 


tion on the set 


K (x) } • 


This is, of course ^ the simplification used 


in a P^-approximation and would not yield exact moments. The familiar 
case of M = 0 poses no such difficulties and one can readily find 
the exact moments. 

In considering the collision density moments, let us assume 

that A(n) is an even M degree polynomial in p. It follows that 

Fj^ * 0 for n + j odd. It is also possible to show, from Eq. (2.5*^) 

and the properties of |a j , that = 0 for n > J(M + 1) and 

therefore, with finite M, we can calculate any collision density 

♦ 

moment exactly. For example, we find for the case of M = 2, 

„ „ v F 0 _2 18 A 2 
< x > = o - 2L P + 175 (1 - c)(l - Cf 3 ) 
r o > 

where ^x 2 ^ is the mean value of x 2 and L^, Is the diffusion length 
as calculated by a collision density based Pj -approximation. We 
note that the result of Eq. (2-5*5) is unlike that found in the 
angle -independent mean free path case. In that case M = 0 and we 
find that the second spatial moment is given correctly by diffusion 
theory, i.e./^x 2 ^ * 2L 2 both in the exact calculation and in diffusion 
theory. 
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In passing, we note that the set can be found from 

the set |f^| via the easily derived relation 


♦J- 




( 2 . 5 . 6 ) 


The result of Eq. (2.5*6) does not contradict our earlier assertion 
regarding the problem of finding the flux moments. When the determi- 
nation of the Is approached directly by Eq.. (2.5.2), it is the 

total cross section Which is considered to be an M degree polynomial 
in 4 . When is found using |f^| as in Eq. (2.5*6), the mean free 

path is assumed to be an M degree polynomial in 4 . 


2.6 The Case of Isotropic Scattering 

I 

The stationary state form of Eq. (2.2.1) with the 
additional assumption of isotropic scattering gives the collision 

l 

density equation 


a 


U Mu) F(x,4) + F(x,4) 




♦ 1 

F(x, 4') dp' 


+ S(x, 4 ) 


( 2 . 6 . 1 ) 


In this case of isotropic scattering an angle variable change is sug- 
gested. Specifically, let us define the angle variable u = 4 A( 4 )/A(l), 
measure x in units of A(l), and change dependent variable to F(x,u) = 
g(u) F(x,u) with g(u) = Idp/duj. With the source density change 
S(x,u) = g(u) S(x,p), Eq. (2.6.1) takes the form 

u 5 ^-F(x» u ) + F(x,u) = | g(u) j F(x,u' ) du* + S(x,u) (2.6.2) 

*'-1 
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We shall consider various aspects of the solution of Eq. (2.6.2). 


Iegendre Polynomial Expansion 


following the procedure used in the derivation of Eq. (2.1.10), 
we obtain the coupled differential equation form of the transport 
equation i 


dF 

n n- 
2n + 1 dx 


i dF 

n ^ 1 n+i 
2n + 1 dx 


F n< X > " 2iT 


c 


+ 1 


M*) 




where we have used the Legendre polynomial expansions 


F(x,u) = ^V~^ F n (x) P n (u) 


S(x,u) = 2 g2 2 Ll S n (x) P n (u) 


g(«) = X) g n P n (u) 
n 


It should be noted that the requirement of a symmetric k(p) imposes 
the condition that =0 for n odd. 

The idea of a P -approximation is equally well -applied here. 

a 

For example, the Pj -approximation (diffusion theory) takes the 


(2.6.3) 

(2.6. 4a) 
(2.6.4b) 
(2.6. 4c) 
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f 


usual form (cf., Eq. (2-3*7)) 


1 * Iy o 

J — Ml -eg,) 


dx 


F, 


ifo 

3 dx 


(2.6,5a) 

(2.6.5b) 


where only isotropic sources are allowed. 

Moment Decomposition 

With tbe collision density and source moments defined in 
the usual manner (i.e., by Eq. (2.5.1)), Eq. (2.6.3) can be trans- 
f ormed to the algebraic set 


n 


gj , C g n 

n 2n + 1 


A 


2n 


tt [" < n - £!] 


( 2 . 6 . 6 ) 


We have assumed that A(u) is symmetric which implies that g(u) i6 

symmetric. We also find that rj = 0 for j + n odd, and, for oad n, 

depends only on j and F^.j . Furthermore, we have the interest- 
n n-l n+i 

ing property that the spatial moment Fjj depends only on the set of 
moments |f* , n + i < . Therefore, the calculations of a low- 

order spatial moment requires the specification of a small number 
of the r and the prior determination of a small number of other 


moments. 


-2k- 


i 






As an example , let us calculate, by these methods, the second 

spatial moment of the neutron distribution resulting from a unit, 

plane, isotropic source (at x = 0}. In this case = 5 „ 6. . 

n no jo 

0 0 1 

The moments Fq , F 2 , and Fj are easily determined and are the only 

2 


values required in the calculation of F 0 . We fipa, 
malized second spatial moment (cf . , Eq. (2-5*5))> 


<x2> 



2 1 + 2c g 2 /5 

3 1 - c g 0 


for the nor- 


Normal Mode Expansion 

We shall apply the recently developed normal mode technique 
(C) to the problem of determining the exact and asymptotic solution 
of Eq. (2.6.2). In so doing we shall arrive at an interesting 
mathematical problem the details of which we consider in Appendix B. 
We consider the homogeneous form of Eq. (2.6.2), i.e., S = 0. Trans- 
lational invariance suggests the "ansatz" 


F(x,u) = 4> (v,u) exp(-x/v) 


where we allow the separation variable, v, to be complex. We obtain 
the integral equation 


(v - u) o(v,u) = | v g( 


♦1 

a) l ' 


*(v,u') du' 


We adopt the normalization 

♦1 



♦(v,u) du = 


1 


(2.6.7) 


( 2 . 6 . 8 ) 


(2.6.9) 


( 2 . 6 . 10 ) 
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» 

If we allow solutions of Eq. (2.6.9) to be distributions (in the 
sense of L. Schwartz (S)), then we have 

t 

♦(v,u) = £ V v ^ U ^ + A(v) 6(v - u) 

We assume that g(-) satisfies a Holder condition (M, p. 11) on 
the interval of the real line u e (-1, +l). Any singular integrals 
which might appear are then of the Cauchy type and we define their 
evaluation as the Cauchy principal value (M, p. 26). 

Hie normalization required of $(v,u), i.e., Eq. (2.6.10), 
leads to a specification of allowed discrete values of v in the 
region of the v-complex -plane excluding the line (-1, +l), and 
to a specification of the function A(v) for v e (-1, +l). To aid 
in the analysis of these results, let us define the Cauchy integral, 
G(v), by 

♦1 

G(v) = f da 

2x i £ u - v 

With v^(-l, +1), we find that Eq. (2.6.10) gives 
1 + i n c v G(v) = 0 

which has a set of roots which are distinct. With v £ (-1, +1), 

Eq. (2.6.10) yields an explicit formula for the function A(v) and no 
restrictions are placed on allowed values of v. We find 


( 2 . 6 . 11 ) 


( 2 . 6 . 12 ) 


(2.6.13) 


A(v) = 1 + i n c v G(v) 


(2.6.14) 
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Thus, if we extend the definition of A(v), bb expressed in Eq. (2.6.14), 

t 

to the entire v-plane, we find that the zeroes of A(v) determine 

the set of allowed distinct values of v. Since g(u) is symmetric, 

G (-v) * -G(v), whence, A(v) is an even function of v. The zeroes 

of A(v), therfore, appear in pairs which we label ±y . 

J 

We have found a set of functions of the angle variable u 

indexed by v, |*(v, u)| . There is a discrete indexed set with 

v£(-l, +1) and members characterized by 

P v g(u) 

♦ (±v.,u) = ^ 


2 v T u 
J 


(2.6.15) 


and, a continuous indexed set with vC (-1, +1) and of form given 
by Eq. (2.6.11). The function A(v), which appears in $(v,u) for 
v £(-l,+l), is given by Eq. (2.6.14). Furthermore, the zeroes of 
A( v ) for v £(- 1, +1) establish the set of discrete indices J±Vj| . 

If we assume that g(u) £ 0 for u £ (-1, +l), we may write 
Eq. (2.6.9), with the normalization of Eq. (2.6.10), in the form 

[1 _ h | = £ 

L vj 2 

Let us multiply Eq. (2.6.16) with index v by ♦(v',u) and subtract 
the result from Eq. (2.6.16) with index v' multiplied by o(v,u). 
Employing Eq. (2.6.10) and Integrating over u £ (-1, +1) we obtain 


♦1 

[v - ^]^ £?IT* (v ’ u) * (v ’' u 


) du = 0 


(2.6.16) 


(2.6.17) 
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Hiere is clearly no degeneracy and thus Eq. (2.6.17) may be 
rewritten as the orthogonality relation 

— *{v,u) ♦(v’,u) du = 0 for v # v' 



Bie nature of the orthogonality relation including the case 
V * v' depends (bn whether v is a member of the discrete index set 
or belongs to the continuum. If v is a discrete Index, then 


♦1 


/ irir * (± v u) * (lv i' u) du ■ 6 ji 

S' 

i(*v ) . 4 j / 

J i ( Y U>< 


du 


If v belongs to the continuum, then 

♦1 

f ♦(v,u) ♦(v , ,u) du = B(v - v’) 

-1 

We have found that the set of normal modes, |o(v, u)| , is 
orthogonal, with weight function u/g(u), on the interval u C(-l, +1). 
For the remainder of this section we shall assume that the normal 
modes are also complete in the space of functions which satisfy a 
Holder condition on the interval u £ (-1, +1). In Appendix B we 
shall, in measure, substantiate this hypothesis by demonstrating 
the existence of the modal expansion coefficients. In so doing, we 
shall generalize the interval of completeness to all physically 


(2.6.18) 

( 2 . 6 . 19 a) 

( 2 . 6 . 19 b) 

( 2 . 6 . 20 ) 


relevant canes. 
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Assuming that the normal inodes form a complete Bet on the 
interval u £ (-1, +l), ve have the general solution of Eq. (2.6.2) 
in the form 

F(x,u) = 53 a(v) #(v,u) exp(-x/v) 
v 

* 

where the summation indicates integration over continuous spectra 
when applicable. In many problems we find boundary conditions 
which can be formulated as 

F(o,u) = «• ( u) = 53 a ( v ) *(v,u) for u £ (-1, +l) 
v 


Mid we can use the orthogonality relations to determine the expansion 
coefficients, a(v). In detail, Eq,. (2.6.22) is rewritten as 


♦(u) = £ a ( +v 1 ) + 53 a(-v ) ♦(-v.,u) 


♦ 1 


+ / a(v) $(v,u) dv for u £ (-1, +l) 

-1 

Direct application of the discrete index orthogonality relation, 

Eq. (2.6.19), yields the discrete indexed expansion coefficients, 

♦ 1 


a(±v j) ‘ X 


%t/ iW * (u 


) »(±Vj,u) du 


-1 


Using Eq. (2.6.18) we obtain, from Eq. (2.6.23), 
♦ 1 ♦1 


/- 1 ♦ 1 

J ♦(») du = / du g(u) *^ V ' U -/ a(v' )0(v',u)dv’ 


( 2 . 6 . 21 ) 


( 2 . 6 . 22 ) 


(2.6.23) 


(2.6.24) 


(2.6.25) 
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There appears a doubly Cauchy singular integral and thus the order 
of integration in Eq. (2.6.25) nay not be reversed without due 
caution. The doubly singular term appears as 



We assume that a(v) satisfies a Holder condition for v £(-l, +1) 
and follow the dictates of the Poincare'-Bertrana formula for invert- 


ing the integration order (M, p. 57). We find 




g(v) a (v) 



v* a(v') r uj 

v ' “ u J . i v 


Using Eqs. (2.6.20) and (2.6.26) we obtain the more useful "ortho- 
gonality relation, " 


£ du gTIT ( a(v') ♦(v , ,u) dv’ = t(v) a(v) for v £(-l, +1) 


•/-l 


I(v) = V g(v) 


im) + 


Now, applying Eq.. (2.6.27) to the problem of finding the continuum 
expansion coefficients in Eq. (2.6.23), we have 


a (v) 


1 

iTvT 


£ 


,♦1 

u 


$(u) $(v,u) du 


Jil d. u 
• u 

( 2 . 6 . 26 ) 

( 2 . 6 . 27 a) 

(2.6.27b) 


for v C (-1, +1) 


(2.6.28) 
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2. 7 Pie Green's Function for the Case of Isotropic Scattering . 

As a specific example of the use of the relations just 
developed, let us consider the problem of finding the infinite 
medium Green's function for isotropic, plane sources. In this 
case, the source density, S(x,u), of Eq. (2.6.2) represents a unit, 
plane, isotropic emission of neutrons at a position which we choose 
to designate x = 0, i.e., S(x,u) = g(u) 6(x)/2. Integration of 
Eq. (2.6.2) over a vanishing interval about x = 0 yields the bound- 
ary condition 

u £f(0 ,u) — F(0 ,u)] = g(u)/2 for u C (-1, +1) i (2.7.1) 

ft 

We impose the additional condition that as |xj— * 00 , F(x, u)— * 0 
and express the solution in the form 


F(x,u) «= a ( + v ) *(+v ,u) exp(-x/v ) 
j j J J 




■ I a ( v) *(v,u) exp(-x/v) dv for x > 0 (2.7.2a) 


F(x,u) = - ^ a(-v^) o(-Vj,u) exp(x/vp 

ij 


- / a(v) 4>(v,u) exp(-x/v) dv for x < 0 (2.7.2b) 


The source condition, Eq. (2.7-1), then takes the form of the general 
boundary condition, Eq. (2.6.23). Specifically, we have 


♦ 1 


sfer ^ a < +v .)>* (+v j’ u > + ifer ^ a( - v j ) • ( - v J ’ u) + ifer/, a ( v )*( v -u)dv 


_1 

2 


( 2 . 7 . 3 ) 



| 
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Whence, employing the normalization expressed by Eq. (2.6.10), 


a(±Vj) 


«(v) = 2iTvT f0r v £ +1 ) 


and we have completed the solution of the Green's function. 

let us examine some aspects of the Green's function. For 
simplicity, we assume that there is only one pair of zeroes of A(v), 
±v fr . We shall develop a sufficient condition for this property 
in Appendix A. 3he Green's function is then given by 


F(x,u) » 


*( + v 0 , u) 


r * 1 

«p(-x/v 0 ) + / l [^c xp(-x/v 

Jo 


) dv for x > 0 


*(-v # ,u) 

F(x,u) = - - ^ fr zvr ' i ex P( x /v ( 




’•> -f M 


exp(-x/v) dv for x < 0 


With the definition 


r * 1 

♦ n ( v ) = J *(v,u) P Q {u) du 

and the easily derived symmetry properties* l(-v # ) = - l(+v 0 ) and 
I(-v) = - I(v), we find the collision density moments for the neutron 
distribution from a unit, plane, isotropic source, 

T n v J+1 

*2 = £ [V +v o) + (-l) j ^(-vo)] 

• fl 

+ 2~ jf [*n (v) + ( ~ 1)J f n (“ v )] dv 
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Using Eqs. (2.1«9)> (2.6.4c), (2.6.9) and (2.6.10), we obtain 
a recurrence relation for the set j* n (v)| , 


(2n + 1) v * n (v) » n + (n + l) * n+i (v) + c v 


g, 


n 


and the normalization ♦„(!») = 1. We note that Eq. (2.7.8) 
implies that * Q (v) is an even, or odd, polynomial in v of degree n. 
Therefore, ve have ^(-v) * (-l) n ♦ (+v), and Eq. (2.7*7) reduces 
to 


r> = 


J! 


V c , 

iTnCT ‘n (+v 


r*V. 

•> * / V 

» o 


= 0 


t (v) dv 

if j + n is even 

* 

t 

if J + n is odd 


We have already considered the moments set 


{ F n} ' 


For this 


particular case the I’n} is determined from Eq. (2.6.6) and the 
source condition & jo* 111 P* 883 - 11 ^ we note that the con- 

sistency of Eq. (2.6.6) and (2.7*9) is easily demonstrated via the 
recurrence relation Eq. (2.7.8). Moreover, equating the Fjj and 

moments as derived by the two relations, we find 
♦1 


/ 


l(+vj + / l(v) dv * 1 - c g 


T FvJ 


/ 


v 3 ^ 1 + 2c gj/5 

(T77T7 


(2.7.8) 


(2.7.9a) 

(2.7.9b) 


(2.7-lOa) 


(2.7.10b) 
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From Egs. (2.7«10a) and (2.7-10b) we obtain an explicit expression 
for the discrete index, v v i.e.. 


v 


M 


1 

3 


(i + 2c g^/5) r +1 v s 

“cTeT "io 



(2.7.11) 


For c < 1, v # is real and is interpreted as the exact asymptotic 
diffusion length (here, measured in units of X(l) ). It should be 
noted that the integral terms in Eg. (2.6.39) depend on c and | g n | 
via the dependence of l(v) on these parameters (cf. , Eg. (2.6.2710). 


III. REMARKS REGARDING APPLICATION OF THE THEORY 

We shall develop a limited number of considerations relevant 
to the application of the theory presented in Section II. lliese 
remarks are intended as a brief illustration of possible methods 
of application of the present theory to physical problems. Many 


l 



interesting calculations are possible, and with the accomplishment 
of experimental measurements of neutron distributions in the types 
of media tender discussion, many comparisons of theoretical and 
experimental results would be profitable. 

Ud certainly require methods of determining the proper varia- 
tion of mean free path if these mathematical formulations are to 
be applied to physical problems. In this section we shall discuss 
the general types of heterogeneity toward which the current theory 
applies. We shall detail a simple method, using known diffusion 

lengths, to specify the angular dependence of the mean free path 

* 

for a particular type of heterogeneity. 

3-1 Types of Heterogeneity 

As mentioned earlier, the motivation of the present 

effort is the establishment of a method of homogenization of regular 

arrays of vacuum channels for the purpose of neutron diffusion 

calculations. We also imposed the necessary restriction that, in 

* 

general, the type of heterogeneity considered should yield two 
characteristic orthogonal directions. As an example of the caution 
which must be exercised in application of the theory, let us consider 
a type of heterogeneity which, at first approach, appears to satisfy 
the necessary requirements, but which actually is unsuitable for 
these methods. Specifically, we examine the case of a periodic 


slab array of scatterer and vacuum. This heterogeneity exhibits 
two characteristic orthogonal directions} perpendicular to slab, and 
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the directions in the plane of the slab. Moreover, the direction 
perpendicular to the slabs (transverse to slab "channels") yields 
considerations which are algebraically easily accomplished. If A(p,x) 
represents the mean distance traveled to a collision by a neutron 
located at a position x to the left of the right-hand-face of a 
slab of scatterer, traveling with direction cosine p relative to 
the slab perpendicular direction, we find 

T 

~ exp(-x/X u) 

A(p,x) = A + for x < T , p > 0 (3.1.1) 

“ l-exp(-T s /X s p) 


In Eq. (3. 1.1), A is the mean free path in the scatterer material 
s 

which has slab thickness T , and T is the vacuum slab thickness. 

s v 

4 

For the homogenized medium we require a function A(p) which, it 

would seem, should be a "suitable" average of A(p,x). For the case 

of isotropic scattering, the average 

T 

s 

y n (x) dx 

o i 

4 ) “ «ji 




* 0 (x) dx 


(3.1.2) 


is clearly indicated. In Eq. (3-1.2), t Q (x) represents the actual 
angular integrated neutron flux. We note that, in the present case, 
ty Q (x) can be found. Far from neutron sources we have ^ 0 (x) — ► exp(x/L g ) 

where L g is the "Itsymptotic" diffusion length in the scatterer material. 


t 
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A note in passing: If T /L « 1, then ijr (x) is approximately 

*" B 6 O 

constant and Eq. (3- 1.2) gives the result 

Mu) - \ (1 + T^Tg) (3.1.3) 

which is the "simply homogenized" parameter. Of course, the condition 

T /L « 1 should yield the homogeneous limit. 

If x now represents the direction perpendicular to the slabs 

we have the asymptotic result ^ (x) — ►exp(-x/L ) when the position 

o s 

x falls in a scatterer slab and * (x) is a constant when x falls 
in a vacuum slab. Ifce "best fit" to this flux, for the homogenized 
medium, is t 0 (x) — ► exp(-x/L) where L is the simply homogenized dif- 
fusion length, i.e., L=L (1 + T/T). This result would be obtained 

S v 6 

if Eq. (3*1.3) were used. In this particular case, we have the 
situation that a calculation based on an angular -dependent mean free 
path yields results that are less representative than the simply 
homogenized calculation. It is expected that, in the orthogonal 
characteristic direction (i.e., in the plane of the slab), use of 
an angular-dependent mean free path is indicated. 

The case of a calculation in the slab perpendicular direction 
for a periodic slab array is certainly excluded from the present con- 
siderations. Moreover, one should feel no motivation toward develop- 
ing a theory for that case since it is easily treated by a standard 
method, i.e., change of position variable to "optical thickness." 


We shall new consider the details of a macroscopic-parameter-based 
calculation for a heterogeneity for which the present methods were 
clearly intended, i.e. , a regular array of cylindrical vacuum channels. 

3*2 Cylindrical Channels in a Regular Array 

Let us consider a regular array of vacuum channels of 

cylindrical cross section. With every vacuum channel of cross 

sectional area A^. we associate a cross sectional area of scatterer 

material A such that V = A Jk is the ratio of vacuum volume to 
s vs 

scatterer volume characteristic of the medium. We shall label the 
axial, or longitudinal, direction with x and direction cosine 
and the radial, or transverse direction with y and direction cosine 
ij. Due to streaming along channels we expect different diffusion 
properties in the x and y-directions and both of these cases to be 
different than the simply homogenized diffusion. The simply homo- 
genized mean free path is given by 

\ = \ (1 + V) (3-2.1) 

7 

Before presenting a specific method for obtaining a representative 
A(u), let us make some general comments regarding the features of 
such a calculation. It is clear that we have chosen to consider 
only two representative directions) axial, or x-direction, and 
transverse, or y direction. It is also clear that in the present 
lattice the actual description of a straight line path in the 


» 



-38- 




Uv . . 


r 


transverse direction starting from a point in the scattering material 
depends upon both the azimuthal angle about the x-direction and 

the particular position in the scattering material relative to, 

* 

say, the center of a vacuum channel. A "suitable" averaging tech- 

niqua must be employed. Furthermore, ve encounter the same difficulty 

when considering a description of an axially-oriented path. Let 

A (u) and A (tj) represent the angular -dependent mean free paths 
x y 

with respect to x-direction diffusion and y-direction diffusion. 

The following constraints on the "suitable" averaging technique seem 
intuitively reasonable* 


(i) 


(ID 


Hie average mean free path based on axial and 
transverse directions should be equal, i.e., 

hi +1 > 


f \(u) dh = f 

-i •'-i 


* y (n) dT \ 


Hie axial mean free path in the transverse direction 
(i.e., at 4 = 0 ) should be equal to the transverse 
mean free path in the transverse direction (i.e., at 
q = i l), i.e. , 


yo) . y± 1) 


(iii) Both axial and transverse mean free paths should be 
symmetric, i.e., 


* x (h) = \(“h) 


* y (n) * \,(-n) 


(3-3.2) 


(3-2.3) 


(3.2.4a) 


(3.2.4b) 
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A possible method of obtaining k(p) is to find, as a function 

of starting position in the scatterer, the mean free path length 

« 

traveled in all directions- Then, upon "suitably ' 1 weighting this 
quantity (according to whether A^(p) or 7 y(q) I s desired) an average 
yields the angular-dependent mean free path- In even the simplest 
lattice this is a geometric task of considerable magnitude. Here, 
for the sake of brevity, we shall take an alternate, albeit certainly 
less self-contained, route. We shall assume that we have given 
certain macroscopic diffusion parameters, such as diffusion length, 
and use the general constraints of Eqs. (3-2.2), ( 3 . 2 . 3 ) an d (3*2.4) 
to obtain a representation of the mean free path which yields the 
given parameters. To be specific, let us assume that A^(p) and 
>y(Tj) are even quadratics of the respective variables. Thus, in 
terms of the Legendre polynomial expansion 


A 

x 


(w) 


A + A P (p) 

XO X2 2 


(3.2.5a) 


A (tj) = A + A P (p) 
y 1 yo y2 2 


From Eq. (3.2.2), we obtain A „ = A , and this result used in 

xo yo 

Eq. (3*2.3) yields "h =-7^/2. Therefore, in terms of the two 
unknowns, A$ and >2 , Eq. (3-2.5) may be reformulated as 


A (p) = A + X P (p) 

x v 0 2 2 y ' 

\(i) - P 2 (n) 


(3-2-5b) 


( 3 . 2 . 6 a) 


( 3 . 2 . 6 b) 



The arguments used here with respect to the mean free 
also apply to the determination of the total cross section, 
ve expect the general constraints: 


+1 


+1 


(i) J O x (n) <J y (Tj) dtj 


path 

Thus, 


(ii) o x ( 0) = o y (± 1) 


(iii) a (h) = o (-n) and a (rj) = a (-tj) 

A A Jr Jr 

If the total cross section is assumed to be an even quadratic, 
then in terms of the two unknowns, o 0 and a 2 ve have 

a x (u) = a 0 + a 2 P 2 ( 4 ) 

a y (l) = a 0 - | P 2 (Tj) 


For the remainder of this discussion ve shall assume that the 
neutron collision density is used as dependent variable and thus 
the mean free path Is the relevant parameter. One can equally well 
apply these considerations to the neutron flux and total cross section. 

From Eqs. (2.3«10) and (3*2.6) ve have the results 


(\ + 2* /5) J 

\ ■ 3(1 ' "- OU - or,) <3 - 2 - Ta) 

l2 - yp ; 

y _ 3(l - c)(l - cf! ) 


(3.2.7b) 
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We also have 


- cHl - ci 


( 3 . 2 . 8 ) 


From Eqs. (3.1.10) and (3.1. 11) ve obtain 

. <L> 

X L 

s s 


(3.2.9a) 


< L> = Lj3 + 2L/3 


(3.2.9b) 


^ 1 
K * 




(3.2.9c) 


We can use measured values of L and L , or other theoretical treat- 

x y’ 

meats, to find and Iy in order to determine \/^ B and ^ 2 /^ B • 

For example, if ve use Behren's theoretical formulation (B), 


&■ 


1 + 2V + 


+ 2RV 


- 1 + 2V + 


+ RV 


where R is the ratio of the vacuum channel radius to X . In 

s 

Fig. 3 ve present L. and. L/L g as a function of R > RC(0,5), 
for the cases V *0.5, 1.0 and 2.0 as determined by Eq.. (3.2.10). 
Then, in Fig. 4 ve have, for the same values of R and V, the results 
for X Q /X s and X^/X s based on the curves in Fig. 3* 


-ha- 


lt should be 

is Eq. ( 3 . 2 . 10 ) are 
<x?> /2 and<y2>/2 

<x 2 >* 


noted that what we refer to as I? and I? 

x y 

actually calculated by Behrens (B) as 
and, via Eq. (2.5.5), we have 

2L a 18 ^ 

Z \ 1T5 (1 - c)(l - cf,) 

* N 


< y 2 > 


21? 

y 


18 >> 

i ?5 (1 - c)(l - cf,) 


2 

If Tsg > 1 the validity of Fig. k as a relevant representation for 
?v(p) is questionable. However, truncation of X(p) at a quadratic 
would. In that case, also be of questionable usef ullness. 


(3.2.11a) 

(3.2.11b) 


IV. SUMMARY 


We have developed the mathematical formulation of a new 
approach to the homogenization of certain types of heterogeneous 
media (such as a regular array of vacuum channels) for the purpose 
of neutron diffusion calculations. The new method is based on the 
inclusion of an angular -dependent mean free path in the theory of 
neutron transport. In the present effort, calculations are restricted 
to media with plane symmetry and monoenergetic neutron theory is 
employed. Extension to energy -dependent theory and to other sym- 
metries would probably follow the general lines for the familiar, 
angular-independent case without significant additional complication. 
However, it seems clear that the requirement of the existence of two 


I 


orthogonal characteristic directions in the development of the 
angular depeifLence of the mean free path must be imposed. 

We have found that a neutron flux based theory and a collision 
density based theory can lead to significantly different results 
when lew-order approximations, such as diffusion theory, are employed 
in the solution of the transport equation. For the case of isotropic 
scattering, the normal mode technique is applicable, and exact, closed- 
form solutions can be determined. 

Evaluating the results implied by the present theory with 
respect to measurements is impossible. There is a current lack of 
pertinent experimental results for neutron distribution description 
in the relevant type of media. 
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APPEND IX A 


THE FUNCTION A(v) 


We have found previously that the zeroes A(v) for 
v £ (-1, +1) are the discrete set of normal mode indices and that 
they appear in pairs, ±v j- Let us now discuss the number of these 
allowed discrete indices. To this end, and for relations which sire 
useful in Section 2.J, we turn to a brief study of the general 
properties of the function A(v) as defined by Eqs. (2.6.12) and 
(2.6.1V), i.e.. 


A(v) =1+1 n c v G(v) 

+1 


G(v) « 


_JL f e( u ) 

2« i / u - v 

-1 


du 


In terms of the set > as defined in Eq,. (2.§. Vc), A(v) may 

be rewritten 

A(v) = 1 - c v £ g Q Q n (v) 


( A. l) 



A-2 


where Q (v) is a Legendre function of the second kind defined for 
n j 

the entire v -plane by an extension of the Neumann formula (H, p. 51) 
to include v C (-1> +l}> i*e., 

1 f 1 P n (u) 

KM-il rrn au 

-1 

with singular integrals evaluated as the Cauchy principal value. 

For large v, v Q 0 (v) varies as v~ n . Thus, A(v) is bounded for 
large v. Furthermore, the Q^v) are analytic in the v-plane exclud- 
ing v C (-1, + l) and, therefore, A(v) is analytic in this same region. 
We use the contouy illustrated in Fig. 2 and the argument theorem 
(T, p. Il6) to establish the number of zeroes of A(v) in the region 
(-1, +1). Since the zeroes of A(v) appear in pairs, we denote 
the number of zeroes by 2J. The argument theorem applied here yields 

4n J = change in arg A + (u) on C + change in arg A~(u) on C 

We have assumed that g(u) satisfies a Holder condition on 
u C (-1, +1) and therefore G(v) is a Cauchy integral. We apply the 
Plemelj formulae (M, p. 43) to find the limit values G ± (u). We find 

G ± (u) = G(u) ± g(u) 

where G + (u) and G”(u) refer to the limit values of G(v) as v approaches 
u from above and below the real line respectively. From Eq,. (A. 4) 
we obtain the limit values 

A*(u) = A(u) ± c u g(u) 


(A.2) 


(A. 3) 


(A.4) 


(A.5) 



A-3 




New, A(u) with uC (-1, +1) is a real function (with singularities 
at u = ±1), and, we have A(o) = 1 and g(u) is a symmetric function. 
Whence, we obtain the relations 

arg A (u) = — arg A~(u) 
arg A + (o) *= arg A~(o) = 0 
A + (u) = A~ ( -u) 

Hiese results used in Eq. (A-3) yield the number of pairs of 
zeroes, J, in terms of the single angle arg A + (+l), i.e., 

j = i arg a + (+l) 

It should be noted that Eq. (A.3) contains the implicit require- 
ment that A + (u) = 0 for uC (-1, +l). This assumption is not com- 
pletely necessary, however, it probably applies to most cases of 
physical interest and its application greatly simplifies these con- 
siderations. 

We shall develop a sufficient condition for J = 1 in the case 
that g(u) is an N degree polynomial in u, i.e., 

N 

«(“) - £ 6 n P„(u) 

n=o 

We note that the Legendre functions, Q n (v), can be expressed as 
Q n (v) = P n (v) Q 0 (v) - Wq.jCv) 

* Q q (v) = arc tanil v ’ v ( -1 > +1 ) 

- = arc tanh i, v ^ (-1, +i) 


(A. 6a) 
(A.6b) 
(A. 6c) 


(A.7) 


(A. 8) 
(A. 9a) 


(A. 9b) 
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where W (v) is an even, or odd, polynomial in v of degree 
n-i' 

n - 1 (H, p. 51). In these terms A(v) is rewritten as 


A(v) * 1 — C v Q 0 (v) 


E 


N 

P n< v) * E ®a w n-.< v) 
, n*=o 


We also have as u — *+ 1, Q 0 (u)— *♦ od and P n (+l) = 1. Clearly, 

W„(+l) is bounded, and thus, if 

N 



•n=o 


then as u — + 1, A(u) -» - 00 . From Eq. (A.5), In the present 
case, we have * 


N 

A + (u) = A(u^+ u 6 n P n (u) 

n=o 


Therefore, we may conclude the following: If Eq. (A.ll) holds 

and. In the range u G (0, + l), 

N 

2 8 » P n< u) > 0 
n=o 

then arg A + (+l) = * and we have the desired result, J = 1. Let us 
stress that Eqs. (A.ll) and (A.13) give a sufficient, not necessary, 
condition for the number of pairs of discrete indexed normal modes 


(A. 10) 


(A.ll) 


(A. 12) 


(A. 13) 


to be unity. 


APPENDIX B 


A REIZVANT HILBERT FROBIZM 


In Section 2.6 the existence of the modal expansion 
coefficients |a(±Vj), j = 1, 2,***, Jj a(v), v £(-l, +1)| was 
assumed. Moreover, the orthogonality relations are based on the 
whole angle range u £ (-1, +l) and thus only provide a means of 
determining expansion coefficients for the case of a boundary 
condition given over all angles. By reducing the problem of 
finding expansion coefficients to the solution of an inhomogeneous 
Hilbert problem, we find that one can demonstrate the existence 
of expansion coefficients, and, prescribe a method for determining 
the value of the coefficients for problems involving all physically 
relevant boundary conditions. We follow closely the techniques 
elegantly described by Muskhelishvili (M). 

Reduction of Transport Problem to an Inhomogeneous Hilbert Problem 


We shall, in general, encounter transport problem boundary 

conditions of the form 

J J 

♦ (u) = ^ ft ( +V j) *( + v y u ) + ^ ^ a(-Vyu) *(-Vj,u) 


J'l 


J=1 


B 


+ / a(v) ${v,u) dv for u£(a,B) 

“a 


(B.l) 


where -1^CI<P< + 1. Let us suppose that by some method we 
are able to determine the set of discrete indexed coefficients 
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|a(±Vj), J = 1, 2, •••> and define 


J J 

*'(u) = *(u) - a(+Vj) *( +v j> u ) “ a (-Vj) 

J=1 J=1 

We then have an integral equation for a(v), vC( 3 >P)> i-e.. 


f. 


a(v) *(v,u) dv - $ ' ( u) for u.C(2> p) 


Using the derived form of «(v,u) Eq. (2.6.11), we obtain 


A(u) a(u) + |g(u) f ^ dv = 4>'(u) for u C(3, p) 

*4 


From Eq. (A. 5 ) 


U = 2n~T [ A+ ^ u ^ “ A "(u)] 
A(u) = | jV(u) + A’(u)J 


and therefore Eq. (B.4) may be rewritten as 

1 rA 


g- [A (u) + A“(u)jua(u) + ^A + (u) - A~(u)J A(u) 


t 


= u t ' (u) 


A(u) = 


2k i 


6 

/* v a(v ) 

I V - u 


dv 


(B.2) 

(B.3) 

(B.4) 

(B.5a) 

(B.5b) 

(B.6a) 


for uC(«, P) 


(B.6b) 




B-3 

We have assurance that A(u), as defined in Eq. (B.6b), exists if 
a(u) satisfies a Holder condition on uG(3, 3). For the moment, 
let us assume that this condition is fulfilled and define the 
Cauchy Integral, A(v), over the entire v-plane, 

fl 

A(v) = r-i— r f - U - - a -^- U - ^ - du 

v ' 2n i / u - v 

Ihe PlemelJ formulae yield the limit relations on the line u C (3, 3) 
A + (u) — A (u) = u a(u) 

A + (u) + A (u) = 2A(u) 

The results of Eqs. (B.5) and (B.8) applied to Eq. (B.6) 
give the alternate form 

A + (u) A + (u) — A (u) A (u) = u <>'(u) for u C (a, 3) 

We have assumed that A~(u) ^ 0 for uC(2, 3). With this 

restriction we can easily transform Eq. (B.9) to the form of a 

boundary condition for an inhomogeneous Hilbert problem on an arc 

\ 

(M, ch. 10). Restating the problem of determining a(u), in these 
terms: Find the sectionally analytic function, A(v), vanishing at 

infinity, with boundary condition on the line uC(ct, 3), 

A + (u) = a» + 

A (u) A (u) 


(B.7) 

(B. 8a) 
(B.8b) 

(B.9) 


(B.10) 
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We note that the assumptions on g(u) and A~(u) imply that 
A~(u)/A + (u) is a function satisfying a Holder condition and not 
vanishing on u £({X, 3), and, if we assume that the angle boundary 
condition, *(u), satisfies a Holder condition and a(±v^) exist, 
then u ♦' (u)/ A + (u) satisfies a Holder condition on u C (&, 3). 

Let us help clarify our procedure by summarizing. If we 
assume (what we wish to prove) that &(u) satisfies a Holder condition 
on u £(<*,£), then the integral A(v), defined by Eq. (B.T), is 
of the Cauchy type. Now, Cauchy integrals are sectionally analytic 
functions witli boundary the line of integration. Specifically, if 
(3,3) is the line of integration: 

(i) f A(v) is analytic in v-plane excluding (a, 3). 

(ii) A(v) approaches well-defined limits as uC (a, 3) 

is approached from either side of (a, 3) with possible 
exception of the end points, a and/or 3 . 

(iii) Near the end points, A(v) satisfies the conditions 

|A( v) j < as v— ►a 

v - a| 


|A(v)| < 7h aB v— *-3 

|v-3| 

where a, b, A and B are real constants, and a < 1 


and b < 1. 
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Moreover, A(v) vanishes as jv| — * qd. We have transformed the integral 
equation for a(u) into the boundary condition Eq. (B.IO) which is 
the form of an inhomogeneous Kilbert problem boundary condition. 

Thus, we have reduced the original transport problem to an inhomogeneous 
Hilbert problem. If we can find a solution, A(v), wh^ch introduces 
no physical ambiguity, then our assumption of the existence of a(u), 
uC (3, 0), will be substantiated. 


Solution of the Hilbert Problem 

In terms of 6(u) * arg A + (u), we have A~(il)/A + (u) = exp(-2i 9(d)) 
and the Hilbert problem boundary condition (cf. , Eq. (B. 10)) 


A + (tt) - exp(-2i 9(u) ) A - (a) + ^ for uC (3, 0) 

A + (u) 


Since A(v) must also vanish as j v j — ► co , the solution is (M) 


A(v) 


H(v) 
2n i 


/ 


u »* (u) 


(u - v) A + (u) H + (u) 


du 


(B.ll) 


(B.12) 


where H(v) is the fundamental solution of the associated homogeneous 
Hilbert problem and is given by 


H(v) = (3 - vr e(3) /* o - v) e(3)/n e 6(v) 
The Cauchy integral ©(v) is defined by 



(B.13) 


(B.lM 


1 
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Providing k = 0(£)/n - &{%)/* is a positive integer, ve have the 
k additional requirements 



(B.15) 


These additional requirements are a necessary feature of the solution. 

It should be recalled that the function $'( 0 ), u C (a, fi), is not 

completely specified, i.e. , the discrete indexed expansion coefficients, 

a(±v.), in Eq. (B. 2) are, as yet, unknown. For the general pro- 
J 

blema considered later, it will be demonstrated that, in each case, 

the k requirements axe necessary and sufficipnt for the complete 

> 

specification of all discrete and continuum expansion coefficients. 


Application of the Hilbert Problem Solution 

Plane symmetry transport problems fall into two general 
categories: 

(i) Infinite media problems with full -angle -range 
boundary conditions (such as the Green's function 
solved in Section 2.7)* 

(ii) Half -space media problems with half -angle -range 
boundary conditions (such as albedo or Milne type 
problems). 

Combinations of the solutions of these type^ problems lead to the 
solution of cases with finite media (slabs). For full-range boundary 
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conditions, the orthogonality of the normal modes provides a direct 
method for determining expansion coefficients. The solution of 
the Hilbert problem in these cases demonstrates the existence of 
the coefficients and thus partially supports the completeness 
hypothesis. For half -range boundary problems, there are no apparent 
orthogonality properties of the normal modes. In these cases, the 
solution of the Hilbert problem not only provides proof of existence, 
but also gives a well-defined prescription for the determination 
of expansion coefficients. We shall now outline the application 
of the Hilbert problem solution to the categories of full-range and 
half -range boundary conditions. 

In the case of an "infinite medium, full-range boundary con- 

f 

dition problem, a source condition is usually given at some position, 
which we choose to designate x = 0. For c < 1, it follows that 
F(x,u) should vanish as jxj — » oo. Ihus, the general form of solu- 
tion is that given in Eq. (2.J.2). The source condition can be 
f ormulated a6 


o(u) « 


£ 

«j-i 


a ( +V j) 


« , ( + Vj,u) 


J 

^ ] a(-Vj) »(-v,,u) 

J“1 


r 

+ I a(v) ♦( v, u) dv for u C (-1, +1) 

( J - 1 

Instead of using the obviously indicated orthogonality properties. 


(B. 16) 


let us consider the coefficient evaluation by the route prescribed 
in the Hilbert problem solution. Note that 3 = -1 and 0 = +1. From 
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Eqa. (A. 6) and (A.7 ) , ve have the results 0(-l) = -Jjt and 

©(+l) = Jit. Therefore, in this case, we find that k = 2J and there 

are 2J requirements of the form of Eq. (2.7.15). Specifically, 


n +1 4 « ( ) 

~T T^- du = 0 for n = 0, 1, • * • , 2J - 1 

A+(n\ 

•• \ — / v W 

Eq. (B. 17) provides a sufficient number of equations to find the 

discrete indexed expansion coefficients, a(±v.), j = 1, 2,***, J. 

«J 

The fundamental solution, H(v), is given by (cf. , Eq. (B. 13 ) 



H(v) = (-1 -v)^ (1 - v) J e*^ v ^ 



Thus, A(v) 1 b determined (by Eq. (B. 12)) and we can find a(u) 
for uC (-1, +1) from the limit relation (cf., Eq. (B. 8a) 


u a(u) = A + (n) - A~(u) for u C (-1, +1) 


Since the problem has been completely and unambiguously solved, it 

is clear that the supposition that a(u) satisfy a Holder condition 

is substantiated and we have demonstrated the existence of the 

( 

1 

expansion coefficients. 

For half -space media we consider two types of problems. An 
"albedo problem" is described by a boundary condition at the medium 
surface (x = 0 with medium occupying x > 0) specified for u G (0, +l) 
and the condition that F(x,u) vanish as x — » co . A "Milne problem" 


(B.17) 


(B. 18a) 
( B. 18b ) 


(B.19) 
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is described by a similar boundary condition at x = 0, but with 

F(x,u) — ► *(-v,u) exp(x/v) with v = Vy j = 1, 2, • - • , J, or 
v £(0, +1), as x — ( oo. We have specified these problems as boundary 
conditions on the half -range u£(0, + l)- With obvious modifications, 
the procedure is easily applied to half -space media occupying x < 0 
and boundary conditions on u £ (-1, 0). With the half -space 
occupying x > 0, the general solution of an albedo problem i6 


F(x,u) = ^ * a(+v^) < K +V j> u ) exp(-x/vj) 

J=1 


r 


♦ I a(v) $(v,u) exp(-x/v) dv for x > 0 


(B.20) 


and for a Milne problem. 


o 

F(x,u) = A 4>(-v,u) exp(x/v) + ^ ^ a(+Vj) *( + Vj,u)exp(-x/vj) 

J^l 


r 


+ I a(v) ${v,u) exp(-x/v) dv for x > 0 


In both cases, the boundary condition at x = 0 can be expressed in 
the form of Eq. (B. l), i.e.. 


*(u) * z a(+Vj) *( +v j 


- -f: 


a(v) *(v,u) dv 


(B.21) 


for u C (0, +l) 


(B. 22) 


I 
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Now, Ot * 0 and 0 = +1 and, from Eqs. (A. 6) and (A. l) , Q{ 0) = 0 
and 9(+l) * J«* Thus, k = J and we have the J requirements 

n+1 *' ( ) 

~ LliiL du = 0 for n = 0, 1, • • • , J - 1 

A (u) H (u) 

Tftese are sufficient to determine the discrete indexed coefficients, 

a(+v ), J = 1, 2, •••, J. Hie fundamental solution takes the form 
J 

H(v) = (1 - v) J exp(fe(v)) 



Hie Hilbert problem solution, A(v), vC(0, +l), and the continuum 
expansion coefficients, a(u), u£(0, +l), are found as in the case 
of a full-range boundary problem. Again, we find substantiation 
for the supposition of the existence of the relevant members of 
ja(v)} . Moreover, we find a prescription for calculating the 
expansion coefficients when the use of orthogonality conditions is 



(B.23) 

(B.24a) 
(B. 2^b) 


impossible. 


